xmin = -20;
xmax = 10;
ymin = -25;
ymax = 15;
zmin = 5;
zmax = 46;

N = 50;
xstep = (xmax - xmin) / N;
ystep = (ymax - ymin) / N;
zstep = (zmax - zmin) / N;

[X, Y, Z] = meshgrid(xmin:xstep:xmax, ymin:ystep:ymax, zmin:zstep:zmax);
V = X;

counter = 0;
for i = 0:N
    for j = 0:N
        for l = 0:N
            pol_val = polynom_i6(X(i + 1, j + 1, l + 1), Y(i + 1, j + 1, l + 1), Z(i + 1, j + 1, l + 1));
            V(i + 1, j + 1, l + 1) = pol_val;
        end
    end
    counter = counter + 1
end

figure;
p = patch(isosurface(X,Y,Z,V,0));
isonormals(X,Y,Z,V,p)
set(p,'FaceColor','red','EdgeColor','none');
daspect([1,1,1])
view(3); axis tight
camlight 
lighting gouraud
hold on;

t0 = 0; tf = 200;
x0 = [3.051522 1.582542 15.62388];
[t, X] = ode45('lorenz', [t0 tf], x0);
plot3(X(:, 1), X(:, 2), X(:, 3), 'g');

xlabel('x');
ylabel('y');
zlabel('z');
